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ABSTRACT 

A continuum version of the virial theorem is derived for a radiating self-gravitating accretion 
disc around a compact object. The central object is point-like, but we can avoid the regulariza- 
tion of its gravitational potential. This is achieved by applying a modified Pohozaev-Rellich 
identity to the gravitational potential of the disk only. The theorem holds for general station- 
ary configurations, including discontinuous flows (shock waves, contact discontinuities). It is 
used to test numerical solutions of a model of self-gravitating radiative accretion discs. The 
presented virial theorem should be useful in the analysis of those (possibly radiating) hy- 
drodynamical systems in astrophysics where the central mass and the mass of the fluid are 
comparable and none of them can be neglected. 
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1 INTRODUCTION 

The standard virial theorem applies to a stationary configuration 
of self-gravitating fluid (a star). A textbook level exposition of 
its Newtonian version can be found in |Tassoul| ( [1978} or |Collins| 
( |1978^ . The general-relativistic formulation is still unsatisfactory, 
although preliminary studies are done in 'Gourgoul hon & Bonaz-| 
[zoia|p994) and |Karkowski & Malec| l |2004 ). In this" paper we go 
back to Newtonian gravity, and formulate the virial theorem for a 
steady radiative accretion disc around a central point mass. 

This result emerged from an investigation of stationary radia- 
tive accretion discs done in Mach & Malec'('2012V where a test of 
numerical solutions was required. The choice of the virial theorem 
was natural — one of important applications of its standard version 
is to test numerical models of rotating stars (see e.g.l Hachisu|1986| l. 

The radiation was included in the virial theorem formulated in 
|Anand|fl965^ , but under too stringent assumptions on the radiation 
stress tensor. In this paper the radiation stress tensor appears in the 
equation expressing the conservation of the momentum, but there 
is no need to specify its form. 

The important feature of our work, absent in the former inves- 
tigations, is the presence of the potential due to a point-like cen- 
tral mass. It can cause difficulties in the standard derivation, since 
meaningless (singular) terms do appear. This point requires an ex- 
planation. The classic Newtonian version of the virial theorem can 
be written as 

r + 2 r d^xl-p\V\^ + 2 f d^x^p = 0, 

JvJ 2 J(j3 2 J(|3 2 

where <b is the gravitational potential, p the mass density, U the 
velocity of the fluid, and p its pressure. Imagine that we try to 
apply this version to a system consisting of a stationary toroid 
of fluid and a point mass located at x = 0. Naively, we 
could write p = pf + M(.(5'''(x), where pf denotes the density of 
the fluid only. The corresponding gravitational potential would be 



O = Of + <I>Kep- Here Of denotes the gravitational potential of the 
fluid, and Okcp = -GM^/\x\ is the Keplerian potential of the point 
mass (G denotes the gravitational constant). The fluid related quan- 
tities satisfy AOf = 4;rGpf, while for the central mass we have 
AOkcp = 47rGM^S'^\x). Clearly AO = 4;rGp, as desired. The inte- 
gral 

I d\xp(S> = I dhpf9f+ I d^xpfOKep (1) 

+ f d^ X M,d^^\x)<^f + f rf^xM,5<''(x)OKep 
Jr3 Jk3 

is, however, meaningless due to the last term on the right-hand side. 
In the following we will show that this term does not appear in the 
proper calculation. 

The virial theorem of this paper is presented in the next sec- 
tion. It is then used in order to test the numerical model of a ra- 
diative, self-gravitating disc of gas accreting onto a compact object 
derived in Mac h & Malec| ([2012}. The details of this model are 
given in Section 3.1, and the virial check is presented in Section 
3.2. 



2 THE THEOREM 

Consider a steady self-gravitating configuration of fluid interacting 
with radiation. It satisfies the Euler equation: 

V • (pfU 0V + p + 'P) = -pfVO, (2) 

which we deliberately write in the conservative form (see e.g. |Cas-| 
I tor 2007 1. Here, as in Section 1, U is the velocity, pf the density and 
p the pressure of the fluid. The gravitational potential is denoted by 
O, and the radiation is described by the radiation pressure tensor: 

•p = i r r k kly dQ, 

C Jo J4;r 
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where ly is the radiation intensity, and k denotes a unit vector ori- 
ented in the direction in which the propagation of radiation is con- 
sidered. As usual, V denotes the radiation frequency, and c is the 
speed of light. Here and in what follows, the divergence of any 
symmetric tensor f is understood as a vector field field with com- 
ponents -T)' = dif'' in Cartesian coordinates. The symbol u®y 
is used to denote the dyadic product of two vector fields u and v. It 
is defined as a tensor with components (u ® v)'^ = u'v' . 

We assume that the gas is orbiting around a point mass M^, and 
that the gas and the point mass are isolated from each other (i.e., the 
mass is located outside of the support of pf). We fix the origin of the 
coordinate system so that it coincides with the position of the mass 
Mc. We will work in Cartesian coordinates in this section; they are 
denoted by x', with x = (x' , x^, x'). 

The gravitational potential O can be expressed as O = Okcp + 
<bf. Here the Keplerian potential <I)kcp is given by Okcp = -GM^I\s.\, 
and the gravitational potential of the fluid satisfies the Poisson 
equation: 

AOf = AnGpf. (3) 

We will now use a modified version of the Pohozaev-Rellich 
identity ( |Pohozaevl 1 1 965 [ |Rellich|[l940| l for the above equation. 
Consider a vector field 

w = ((x • V)cl)f + Of /2) VOf - X |VOf p /2 + AnGpx. 

A straightforward calculation making use of the Poisson equation 
l[3j shows that 

V • w = AnG (pfOf/2 + pfX • VOf -I- x • Vp + 3p) . 

Note that Okcp does not enter the above formula at this stage, and 
there are no aforementioned issues with the singularity at x = 0. It 
follows from the Euler equation (|2| that 

PfX ■ VOf -I- X ■ Vp = -pfX • VOfep - X • (V • (pfU ® U)) 
-X • (V • P). 

Moreover, for the Keplerian potential, we have -x • VOxep = Oxep. 
Thus 

V-w = 4;rGpf((0 + OKep)/2-x-(V-(pfU(g)U))-l-3p 
-x • (V • f)) . 

We will assume further that the fluid is finite, i.e. the support 
of the density pf and the pressure p is compact. The reader can 
consult |Mach & Simon|pOT2} for the necessary assumptions in the 
infinite case. We integrate the above equation over a ball S(0, R) 
of radius R, centered at x = and containing the support of the 
pressure and the density. We then take the limit of /? — > oo. The left 
hand side of the obtained equation can be converted into a surface 
term that vanishes, since the gravitational potential Of falls off like 
l/|x| asymptotically. Moreover, the second term on the right-hand 
side can be integrated by parts. This gives the main result of this 
paper: 

= AnG £^ cPx (pf(0 + OKep)/2 + Pf |Up + 3p - x • (V • f)) 

= AttG (£p„, + 2£k„ + 2£,herm + e) . (4) 

Here Epat = ^3 rf^xpf(0 -1- Okcp)/2 is the total potential energy, 
= ^3 d^xpf|Up/2 denotes the bulk kinetic energy, £thcrm = 

Ju3 (Px3p/2 is the internal thermal energy, and E = - j^^ (fxx ■ 

(V • V) describes the interaction with the radiation. 

The potential Okcp enters the above expression twice: once as 

a part of O, and then explicitly. This corresponds exactly to the 



first three integrals appearing on the right-side of equation l|TJ, as 
anticipated. To see this note that the second and third integrals on 
the right-hand side of equation l[TJ are equal. Indeed, 

r (fxM,S^^\x)^f = \ rf^i;AOKepOf/(4;rG) 



X3 



rf^xOKepAOf/(47rG) 



d xOKcpPf. 



Here the surface terms vanish because of the asymptotic fall-off of 
both gravitational potentials. 

Also note that the term E = - Jd'xx- (V -P) can be integrated 
by parts yielding 



E + surface term, 



where 



E = TiP 



Tdyf /„ 

Jo J4;r 



dQ. 



is the radiation energy density. In general, the surface term appear- 
ing here does not vanish for radiating systems. We will not be using 
this form in the following sections. 

The derivatives appearing in the above formulae can be under- 
stood in a weak (integral) sense (to see this, one has to analyse all 
steps of the derivation of the theorem with respect to the character 
of derivatives). This fact is especially important for hydrodynami- 
cal quantities p, pf, U that can be discontinuous. The gravitational 
potential Of is assumed to be smooth, and its derivatives appear- 
ing in the calculation can be understood in the classical sense. The 
obtained result is thus quite general. Contact discontinuities and 
stationary shock waves are allowed in the underlying solution of 
the Poisson-Euler system of equations j2||3j. An introduction into 
the notion of weak or integral solutions of nonlinear partial differ- 
ential equations with a special emphasis put on conservation laws 
such as (|2} can be found in |Evans| ( |2010| l. 



3 APPLICATION 
3.1 Accretion model 

We will demonstrate that the above theorem allows for an easy 
check of the validity of self-gravitating steady accretion discs with 
radiation obtained numerically in |Mach & MalecH2012^ . 

In the following we sketch the construction of the model; de- 
tails can be found in |Mach & Malec| |20T2|. Consider a quasi- 
stationary configuration of perfect fluid accreting onto the central 
point mass M^. Assume that there exist some mechanism that pro- 
duces radiation with the emissivity per unit mass j, and that the 
radiation interacts with the gas via Thompson scattering. It can be 
shown that the radiative flux 



F= r dv f kly 

Jo Jajt 



dQ. 



Jo JAj! 

satisfies V • F = pfj. Moreover, the momentum exchange rate ap- 
pearing in Q is V-y = -(pf/c/c)F, where k is the scattering opacity. 
The energy conservation equation can be written as 



V • (pfU [h + IUIV2 + O) + f) = 0, 



(5) 



where h denotes the specific enthalpy. 

We will now restrict ourselves to polytropic equations of state 
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p = Kp^^ and simple axisymmetric flows. Let (r, 0, z) denote cylin- 
drical coordinates. Assume the velocity of the fluid in the form 
U = Udr + cod^, with U « roj. Define the mass accretion rate 
function M = —InUrpf. It follows from the continuity equation 

V • (PfU) = 0, 

that M = M(z). It is then shown that F can be expressed as a gradi- 
ent of a potential: (/c/c)F = VT. 

The solution of the model is computed from the following 
equations: the integrated Euler equation ^ 

h + (DKep -Ktf + Oc - = C, (6) 

where C denotes an integration constant, the energy equation jSj 
rewritten as 

Avp = — ld^^> + 2ru? + rojdroj) , (7) 
2ncr ^ ' 

and the Poisson equation for the gravitational potential |3]l. Here, 
as usual, the rotation law oj = oj(r) has to be specified a priori, and 
the centrifugal potential (b^ is defined as 

Oe = - J" dr' r'co^r'). 

The radiation mechanism (the emissivity per unit mass j) is 
not specified at the beginning. Instead, we prescribe the accretion 
rate function M and require stationarity of the resulting configura- 
tion. The complete solution can be then computed form the energy 
and momentum conservation equations, together with the emissiv- 
ity j. In this way we are following the spirit of classic papers on 
astrophysical accretion. It has been pointed out by old masters (see 
e.g. |Pringle|198 H that once a steady accretion occurs, one can elim- 
inate the viscosity from the description. The effective model needs 
only M; no specific information on the production of radiation is 
required. 

As a consistency check, it remains to be verified that the con- 
dition C <c roi is valid in the obtained solution. This requirement 
concurs with the fact that solutions of the model do only exist for 
small luminosities, and thus for small M. 

Quasi-stationarity of the model is understood in the usual 
sense; some matter must be delivered to the system from an outside 
reservoir at a small rate M. The central mass grows in time, but 
this growth should be negligible. 

Equations |3|, ([6j and |7| have a scaling symmetry that can be 
used to convert them into a dimensionless form. Assume that the 
disc spreads up to r = Tom at the equatorial plane, and the max- 
imal density of the gas is Pn,ax- The quantity u = GRl^fp^ax has 
the dimension of the potentials. It can be used to define h = h/u, 
•tKep = ^Kep/u, *f = ff/f, <l>c = fc/", = and cb = 

wrout/ V"- Spatial dimensions and the density can be also rescaled: 
X = x/r„ut and pf = p/pma^. The potential Okcp is then expressed as 
OKep = -M,/|x|, where = MJ(p„,.^^rlJ. Equations js}, (je} and 
l|7j can be now reduced to 

h + OKep + Of + 6c - * = C, (8) 

AOf = 47Tp, 

2ncf ^ ' 

where A denotes the laplacian with respect to the rescaled coordi- 
nates X. 




0.2 0.4 0.6 0.8 1 



Figure 1. The density plot of the radiating accretion disc solution. Due to 
the assumed symmetry, only a section through a meridian upper half-plane 
is shown. The rescaled density p is coded in grey-scale. 



3.2 Numerical tests 

Numerical solutions of the model can be found by extending the 
classic Self-Consistent-Filed (SCF) method that was initially used 
to model equilibrium configurations of rotating fluids (stars) jOs-| 
|triker & Mark|I968| l. The basic idea of the SCF method is to convert 
the Poisson equation for the gravitational potential into the integral 
formula, regularise the Poisson kernel by the expansion in the Leg- 
endre polynomials, and iterate it together with the integrated Euler 
equation of the form ([8]l until some level of convergence is reached. 
In the case of the accretion model presented here, there is one more 
Poisson-like equation to be solved, but the overall scheme can be 
retained. Details of this numerical method can be found iu iMach &I 
|Malec|p012l l. 

As an illustration of the virial theorem of this paper we con- 
sider a numerical solution obtained for the rotation law (b{f) = 
cb/r^'^, where o) = const. We search for a disk solution with = I 
and r^lraM = 10"', where rj,, denotes the inner equatorial radius of 
the disk. We assume the polytropic exponent F = 5/3 and the ac- 
cretion rate function 

kM/c= 10"'exp(-(50zf) 

(note that kM/c is a dimensionless quantity). 

The obtained distribution of p is shown on Fig. 1. The mass of 
the disk equals Mf = j^, (Pxpf x 0.294M;.. Choosing r^ut = 10 pc 
and pn,ax = 10"' 'gem"' we get » 1.48 x IO^'Mq. The mass 
of the disk is then Mf » 4.34 x 10' Mq (solar masses). The total 
luminosity can be defined as L = j dS ■ F, where the integral is 
evaluated over the disk surface, and dS denotes an outward oriented 
surface element. It yields L « 2.06 x IO'^Lq, where Lq is the solar 
luminosity. Such a soluion could serve as a model of an accretion 
disc around an ultramassive galcatic black hole. 

For the same solution we get fikin/li'poti ~ 0.446, iiteim/l^potl ~ 
5.33x10-2, and £/|£-po,| « 3.74x IQ-*. Also (j^, rf'xpf 0/2) /E^ot ~ 

0.559, and (/,,, d'xpf9i,,^/2) /E,,, « 0.441. 

The virial check can be performed by computing e = |iipoi -l- 
2£kiii + 2£therm + ^l/l^^potL accordiug to equation l|4|. In our case 
£ X 10"** for the solution obtained on the numerical grid with the 
resolution of 5000 x 5000 points and the multipole expansion using 
more that 20 Legenrde polynomials. Note that this value is quite a 
few orders of magnitude smaller than each of its constituent terms. 
This confirms the validity of the numerical solution. It is also an in- 
dependent check of the absence of trivial mistakes in the derivation 
of the virial theorem of this paper. 
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4 FINAL REMARKS 

We shall warn the reader that a set of functions satisfying the virial 
theorem does not have to be a solution of the Poisson-Euler system 
of equations l|2||3|. Odrzy wolek| ( |2003) presented an iterative ana- 
lytical scheme for obtaining such solutions, where in each step the 
virial theorem was satisfied, although a convergence to the solution 
was not yet reached. 

Despite this, the virial check still belongs to the very few tools 
that allow one to verify the correctness of numerical solutions of 
the Poisson-Euler system of equations (|2||3](. This paper provides 
such a test for a stationary accretion system of self-gravitating gas. 

The virial theorem presented here can be used in the analy- 
sis of those stationary hydrodynamical systems, where the central 
mass and the mass of the fluid are comparable, and thus none of 
them can be neglected. Models of such systems can be found in 
[Hashimoto, Eriguchi & Minier| ( [T995) or |Mach & MaleclpOTT) . 

Magnetic fields can be included in the analysis in the manner 
similar to the radiation. In the astrophysical context this was done 
by |Chandrasekhar & Fermi] ( |1953[ l. As a result, an additional term 
representing the magnetic energy appears in 
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